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In the field of atom optics, the basis of many experiments is a two level atom coupled to a light 
field. The evolution of this system is governed by a master equation. The irreversible components 
of this master equation describe the spontaneous emission of photons from the atom. For many 
applications, it is necessary to minimize the effect of this irreversible evolution. This can be achieved 
by having a far detuned light field. The drawback of this regime is that making the detuning very 
large makes the time step required to solve the master equation very small, much smaller than the 
time scale of any significant evolution. This makes the problem very numerically intensive. For 
this reason, approximations are used to simulate the master equation which are more numerically 
tractable to solve. This paper analyses four approximations: The standard adiabatic approximation; 
a more sophisticated adiabatic approximation (not used before) ; a secular approximation; and a fully 
quantum dressed-state approximation. The advantages and disadvantages of each are investigated 
with respect to accuracy, complexity and the resources required to simulate. In a parameter regime of 
particular experimental interest, only the sophisticated adiabatic and dressed-state approximations 
agree well with the exact evolution. 

PACS numbers: 32.80.Lg, 03. 75. Be, 42.50.Vk, 02.60.Cb 



I. THE MASTER EQUATION 

Many investigations in the realm of quantum physics 
are based on a single two level atom coupled to an ap- 
proximately resonant light field. One major aspect of 
this work is to investigate the quantum mechanical mo- 
tion of such an atom. The work of Steck et al § IS a 
prime example. They investigate chaos assisted tunnel- 
ing by studying the motion of cold cesium atoms in an 
amplitude modulated standing wave of light. Haycock et 
al Q study cesium atoms in an effort to observe quantum 
coherent dynamics. The work of Hensinger et al inves- 
tigates quantum chaos and quantum tunneling 
Proposals in QED to utilize interactions between atoms 
and their cavity for physical realizations of a quantum 
computer and for communication of quantum states 
also require good descriptions of quantum mechanical 
motion. 

The full quantum system here consists of the light field 
and the atom. The evolution of such a system is governed 
by a Schrodinger equation. In most cases however, the 
evolution of the light field is not of interest and so the 
Schrodinger equation can be reduced to a master equa- 
tion. The most general form of this master equation is 
given by 



p = -- [H,p\+Cp, 



(1.1) 



where /! is a Lindbladian superoperator Q]. 



D.AtkinsQgii.cdu.atj 



•Electronic address: 

t Electronic address: H . \ Visenian au . cdu . au 

■f Electronic address: jrahladCemidway.uchicago.edij 



The form we will use here, specific to a two level atom 
interacting with a light field, is 

H = ^ + hda'^a + - \n{x, t)a^ + n''{x, t)a] , (1.2a) 
2m 2 



Cp^r {BJ [(t]p-A [cj] p) . 



(1.2b) 



Here we are only interested in the motion of the atom in 
the one direction. Thus Px is the component of the atom's 
momentum in the x-direction. The detuning of the light 
field is S. The atomic lowering operator is given by cr = 
\g){e\ where \g) and |e) represent the ground and excited 
states of the atom respectively. J7(a:,i) is the complex, 
and possibly time-dependent Rabi frequency operator for 
the light field, which, from here on, will be represented 
simply as il. k and m are the wavenumber of the incident 
photons and the mass of the atom respectively. Planck's 
constant (h) will be set to one for the rest of the analysis. 
The superoperator B describes random momentum kicks. 
It is defined for any arbitrary operator a as 



Ba = 



W{u)e 



ikux ^ 



'ikux 



(1.3) 



where W{u) is the atomic dipole radiation distribution 
function produced by the electronic transition, reduced 
to one dimension . For motion parallel to the direction 
of propagation of the laser light, it is given by 



W{u)^-{l + u') 



(1.4) 



The superoperators and A are defined for arbitrary 
operators c and a as 



J[c\a = cac\A [c] a = - {c^c, a} , 



(1.5) 
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and define the general form of a Lindbladian |^] super- 
operator. 



The enth'e expression of Eq. ( 1.2b ) describes the ir- 
reversible evolution of the system at rate F. It is this 
part of the master equation that we wish to minimize by 
working in the regime 6 3> fimax 3> F (where ilmax is 
the maximum modulus of the Rabi frequency operator). 
The problem with this regime is that when making d very 
large, the full master equation still has to be solved using 
a timestep smaller than S~^. This makes solving the full 
master equation numerically very difficult. 

There are a number of ways to approximate the master 
equation, four of which are investigated here. These are 

1. The standard adiabatic approximation; 

2. A more sophisticated adiabatic approximation; 

3. A secular approximation; and 

4. A dressed-statc approximation. 

Approximation (1) is a standard approach used by many 
researchers in the field, both as a semi-classical treatment 
1^ and as a fully quantum approximation method ^ . 
Unfortunately, this treatment is not valid in the regime 
of the work of |^ . These experiments were performed in 
the regime oi 6 ^ f^max ^ T but where F is the same 
order as fl'^^^^/S. In their work, approximation (1) was 
used but on closer examination, this approximation was 
seen only to be valid in the regime F ^ Q'^^^^/S. We will 
concentrate on the diSicult (F ^ il^^^^/S) regime in this 
analysis. Approximation (2) is one way to correct the 
standard adiabatic approximation for this regime. Ap- 
proximation (3) was proposed in Ref . |]l0| , and was used 
for the formation of quantum trajectory simulations in 
Ref. j^. Approximation (4) is a fully quantum dressed- 
state treatment including the effect of spontaneous emis- 
sion. Previously, only semi-classical treatments, both 
omitting |1^, and including 12 , spontaneous emis- 
sion have been done. We look at all four approximations, 
examining the validity and complexity, as well as the nu- 
merical accuracy (compared to the full simulation) and 
computational resource requirements, of each. 



II. THE STANDARD ADIABATIC 
APPROXIMATION 

The method of adiabatic approximation applied to the 
full master equation serves to eliminate the internal state 
structure of the atom. One reason for wishing to have no 
internal states is so that the system has a clear clas- 
sical analogue [p^ . In performing this treatment, we 
also remove the Hamiltonian term of order 6, removing 
the requirement that the master equation be solved on a 
timestep at least as small as S~^. This adiabatic elimina- 
tion technique is described in many different text books, 
see pSl p3. The result we obtain was first derived by 



Graham, Schlautmann and ZoUer g] . To achieve this adi- 
abatic elimination however, we follow a procedure similar 
to that of Hensinger et al [|j . 

The density matrix can be written using the internal 
state basis as 



P = Ps3 ® \9){9\+Pge «) |g)(e| 
pU ® \e){g\+pee ® |e)(e|. 



(2.1) 



TBpee ~ \: {^'^ Peg ~ Pge^) ~ K^Pgg, (2.2a) 



where the pgg etc are still operators on the centre-of-mass 
Hilbert space L(2)(R). From Eqs. (O), these obey 



Pgg 

Pge 
Pee 



-^Pge- 

'ICpge, 



2 (^Ve 



2 i^Pg^ 



- Pgg^^) + "iSpge 

(2.2b) 

PegCi^) - JCpee- (2.2c) 



Here the kinetic energy term is represented in the super- 
operator 



K-p = i 



2™ 



(2.3) 



If the standard adiabatic elimination procedure is 
valid, we can represent the system by the density ma- 
trix for the centre of mass (com) alone, 



Pcom = Trint[p] = Pgg + Pee, 



(2.4) 



where Trjnt is the trace over the internal states of the 
atom. In this standard adiabatic approximation we fur- 
ther simplify this by noting that large detuning leads to 
very small excited state populations such that pgg ^ p^e 
and thus pcom is approximately pgg. T hus d enoting 

Pcom 

simply as p, we can replace pgg in Eq. (2.2a) just with p. 

Now we require expressions for pg^ and pee i n ter ms of 
p. This is achieved by noting that from Eq. ( ^.2b ), pge 
comes to equilibrium on a timescale much shorter than 
Pgg, at a rate F/2. Thus we set pg^ = 0. Also, if the 
kinetic energy term is much smaller than 5 then it can 
be ignored. This will be the case 5 {p^^/m). This is 
typically true, and so we get 



;Pgg 



17t - l^tpe 



T~2i5 



(2.5) 



This can be substituted back into the equation for pe 
to give 



= -F 



[nvL^Pee] i5[nn\pe 



2(F2-h452) F(F2-f4,52) 



TjmPg^_^ 

+ r2 + 4<52 



(2.6) 



From Eq. (|2.6D, we see that Pee equilibrates on a 
timescale much faster than pgg and so we also set p^e = 0. 
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This time, the kinetic energy term must be ignored com- 
pared with r rather than S. Allowing this approximation 
gives 



Pee + 



{nn^,Pee} iS [nn\p,e] _ npggn^ 
2 (r2 + 4^2) + r (r2 + 452) ~ r2 + 452 • 



(2.7) 



The first correction term on the left hand side scales as 
(f2/(5)2 which in the regime chosen is negligible compared 
to the leading order term. The second correction term 
however scales as {ft^/S) /T. Had we been working in 
the regime F ^ /S, then this term could also be safely 
ignored compared to the leading term. This condition 
however, is not satisfied in the experiments of |3| nor in 
our chosen regime, leaving this term the same order as 
the leading term. This term however was dropped in Ref. 

on the basis that in a more sophisticated approach 
(Sec. Ill), this term does not appear and the correction 
to the final master equation is small Knowing that 
this is an unjustified approximation, but in the interest 
of comparison to currently used techniques we will 
continue to follow this method as others have done. Thus, 
dropping the second correction term we are left with 



r2 + 4(^2 ■ 



(2.8) 



N ow su bstituting first Eq. (2.5) then Eq. ( p.q ) into 
Eq. (2.2a) and replacing pgg with p gives the final adia- 
batically eliminated master equation: 















P ^ r (^Bj 


25 


P- 


-A 


26_ 


p^ - i 



2m 



(2.9) 



where here we have used the fact that (5 ^ 51 ^ F to 
eliminate some of the higher order terms. This master 
equation is of the Lindblad form. 



III. A MORE SOPHISTICATED ADIABATIC 
APPROXIMATION 

As noted, the standard adiabatic elimination method 
described in section || is not strictly valid in the regime 
of the experiments of Hensinger et a/ F ~ ^l'^^^/5. 
There are a number of ways to try and develop a 
strictly valid version of the adiabatic approximation in 
this regime. One way would be to not drop any terms 
without justification and continue to plough through the 
mathematics. Another method which we believe to be 
neater and just as accurate is proposed here by using a 
slightly more sophisticated method similar to that in the 
appendix of p5[ |. 

The basis of the approach is to move into an interaction 
picture with respect to 



(3.1) 



This approach may seem counter- intuitive to most. Usu- 
ally when moving to an interaction picture, it would be 
with respect to a Hq that is already one of the terms in 
the Hamiltonian. In our case, if we investigate Eq. (^.S|), 
we are moving to an interaction picture with respect to 
the opposite of a term in the effective Hamiltonian and as 
such will actually be adding a term to the Hamiltonian. 
The reason we choose to do this is that the problem term 
we encountered in Sec. was in the Hamiltonian for the 
excited state. The potential seen by the excited state 
of the atom is in fact inverted and so the chosen Hq is 
designed to cancel the excited state potential. 

After moving into this interaction picture, we then per- 
form the adiabatic elimination process, and finally trans- 
form back into the Schrodinger picture. This method 
can give a different result because the approximations 
we make in the interaction picture may not have been 
valid in the Schrodinger picture. 

With the unitary transformation operator UQ{t) 

e" 



the interaction picture density matrix is 
p^ul{t)pUa{t). 



(3.2) 



Using this in Eqs. (l.S) gives an interaction picture mas- 
ter equation equation still of the Lindblad form 



P — —i V,p + TCp, 



(3.3) 



but now with an extra Hamiltonian term such that V is 
given by 



V 



2m 2^ 



(3.4) 



Where px is just the component of the momentum in 
the x-direction, transformed into the interaction picture. 
The Lindbladian superoperator C is unaffected by the 
interaction picture because Rabi frequency operator 
commutes with the position operator, X ^ clS well as with 
the state operators, a and cr^. 

Following the same procedure as the standard adia- 
batic elimination process, the equations for the centre- 
of-mass operators can be extracted: 



Paa 


= TBpee- 


\ {^'^Peg - 


- Pge^) - 




-K.pgg, 






Pge 




- ^ {^'^Pee 






+iSpge " 






Pee 


= -Tpee - 




Peg^^) - 




-ICpee, 







AS 



[nn\pgg] 



(3.5a) 



AS 



(3.5b) 



(3.5c) 



where obviously JC uses the interaction picture momen- 
tum operator p^- 

In this more sophisticated approach, we still take the 
trace over the internal states of the atom, letting pcom = 



4 



Pgg + Pee, but now we do not simplify this further and 
simply let p = pgg + p^^ . This gives a master equation of 
the form 



= T [B - I) Pee- I: [n\ Peg]- l:\nrPge] 



^^[nn\p]-icp. 



(3.6) 



Hence we again need to find expressions for pee and pge in 
terms of p. We can achieve this by noting that, as in the 
earlier treatment, pee and pge equilibrate on a timescale 
much faster than pgg. By setting pge = 0, and again 
ignoring the kinetic energy term, we get 



Pge = 



pQt [nn\p^^] {9.\pee] rpf]t 

' 25 + 25 Ai5^ 

r{n\pee} [nn^,{n^,pee}] 



Ai5^ 



8(53 



(3.7) 



Setting Pee = and ignoring the kinetic energy term 
allows us to solve for pee- In the more sophisticated adi- 
abatic treat ment , instead of getting an equation of the 
form of Eq. (|^, we get 



TZpee = 



npfl'' 



(3.8) 



In Eq. (|3.8| ), the superoperators T and TZ are defined 



by performing the opposite unitary transformation. This 
leaves the final master equation for the more sophisti- 
cated adiabatic elimination treatment as 



= r 









-A 


n' 






25 


P- 







pI 



2m 4(5 



16(53 



(3.13) 



Notice here that by explicitly accounting for the term 
in r2^/r(5, the master equation derived is the same as 
that of the standard adiabatic approach with an extra 
potential term of order Q!^ /5'^ . This term is very small 
and as such the statements made to justify the standard 
approach Q were correct in that the adjustment to the 
final master equation is small. The more sophisticated 
treatment however, although more algebraically intensive 
to produce the initial master equation, requires very lit- 
tle extra effort than the standard adiabatic treatment to 
simulate. More importantly though, the master equation 
derived in the more sophisticated adiabatic approach is 
valid in the regime Q.^^^^/ 5 ^ T. 



n 



J- + higher order terms. 



(3.9) 
(3.10) 



In the regime chosen (F ~ ^1^/5), the second term of 
Eq. (3.9) is of order 1 and so cann ot be ignored compared 
to the leading term. In Eq. ( 3.10 ), the higher order terms 
are of order much smaller than 1 and even smaller than 
ri^/(53, which is the smallest order kept by making Taylor 
approximatio ns t o get the expression for pge- Thus to 
simplify Eq. (|3.8|), we act on the left of both sides with 
the superoperator JF~^. Then to leading order we have 



Pe 



npn'f 

4^2 ■ 



(3.11) 



As expected, the interaction picture chosen produced a 
term in pgg to counteract the unwanted term in pee- No- 
tice here that we come to essentially the same result as in 
the standard adiabatic treatment in Eq. (|]^) but with- 
out making any unjustified approximations. 

No w substituting Eq. (3.7) and Eq. (3.11) into 
Eq. (3.6), after simplification, leads to the final inter- 



action picture master equation 







p- A 








(^BJ 








16(53 



-K-p. 



(3.12) 



Finally, all that remains is to transform the interaction 
picture master equation back to the Schrodinger picture 



IV. A SECULAR APPROXIMATION 

A secular approximation to the full master equation 
is quite different from any adiabatic approximation. A 
secular approximation does not totally remove any de- 
pendence on the internal state of the atom. It only elim- 
inates the coherences between the internal atomic states. 
The secular approximation was also used by Hensinger et 
al alongside the standard adiabatic approximation. 

There are a number of ways to derive a secular approx- 
imation to the full master equation. One such method 
has been performed by Dyrting and Milburn ||l^ but this 
method is quite complicated. A much simpler method 
which produces the same approximate master equation 
is based on the technique that is applied for the standard 
adiabatic approximation. 



We start with the Eqs. ( p^ for p 
abatically eliminate Peg 



etc, and then adi- 
However, instead 



as in Sec 

of also trying to solve for p ee, we just substitute the ex- 
pression for Pge, Eq. (| 2.5p, in to the equations for pgg, 
Eq. ( 2.2a ) anti pee, Eq. (|2.2q ). This eliminates the co- 
herences between the excited and ground states while 
still keeping much of the original master equation. One 
advantage of this is that it allows comparison to an anal- 
ogous 2-state classical model ||l^. More importantly to 
us, this approximation eliminates the evolution at rate 5 
as is necessary to simplify the simulations. We are left 
with the following equations for the ground and excited 
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state density matrices 









-A 








26 


Pee 


26 





4(5 



'Pgg 



ICpgg, 







Pgg A 












Pee^ 



~46 



■,Pe 



JCpe 



(4.1) 



(4.2) 



The final master equ ation is constructed by recombin- 
ing the equations (4.1) and (4.2) to produce an equation 
which reproduces these equations for pgg and pee while 
also only giving rapidly decaying terms for pge and Peg- 
This final master equation for the secular approximation 
is 



2m 



V{BJ[a]p-A[a]p) 











A 






+r 




26 


P- 


26 


^) 










A 






+r 




26 




26 





(4.3) 



where tr^ is just the Pauli spin operator a^a — aa'^ . 

To compare this master equation with those we have 
already seen, the Hamiltonian terms derived here are the 
same as those der ived in the standard adiabatic approx- 
imation, Eq. ( ^.9[ ). The spontaneous em ission term ex- 
actly as in the full master equation, Eq. (1.2b), remains, 
while two extra jump terms involving a state change with 
no spontaneous emission have been derived. 

There are no apparent problems in this derivati on, or 
that injioj. Nevertheless, as we will discuss in Sec VIII , 
Eq. (4.3) doe s not give accurate results in comparison 
with Eq. (p]), in the regime VL^ ^ T6. 



V. A DRESSED-STATE APPROXIMATION 

A semi-classical dressed-state treatment of atomic mo- 
tion was put forward by Dalibard and Cohen-Tannoudji 
. The dressed-state approximation used here is a fully 
quantum version of that treatment. 

The states we have been using for a basis so far, \g) and 
|e), are called bare states. We can also work with another 
basis of position-dependent states which we call dressed- 
states. These dressed-states are derived by consi dering 
the Hamiltonian of the full master equation Eq. ( 1.2a ) 
and ignoring the kinetic energy component to get 



V = 5cf^a 



{VLcF^ + n'^a). (5.1) 
The diagonalization of the Hamiltonian yields the form 
V = Ei{x,t)a^a + E2{x,t)aa\ (5.2) 



where Ei and E2 are the eigenenergies of V. The corre- 
sponding eigenstates are the position dependent dressed- 
states we will be considering here, labeled |1) and |2). In 
terms of these dressed-states, a is given by |2)(1|. The 
energies and states are 



Ei{x,t) 
E2{x,t) 

|1> 
|2) 

where is defined by 



i(<5- 
\{5- 



A), 
A), 



suiO\g) - 
cos 6* I g) 



smf 



cost 



V2 ^ 
A + (5 



V2 



6^ + Q.n'^) ' 



(5.3) 
(5.4) 
(5.5) 
(5.6) 

(5.7a) 
(5.7b) 



and A = (52 + rj^2t) 2 ^ 

Now, we use this basis to get equations for pn, pi2, P21 
and P22- Without including the spontaneous emission or 
the kinetic energy, these equations are simply 



P]k = [Ej[x,t)pjk - p.jkEk{x,t)\ 



(5.8) 



We do however want to include the spontaneous emission 
in our treatment and so we must assess the effect of the 
raising and lowering operators for the bare states (cr = 
\g){e\ and a'^ = \e){g\) acting on the new basis states (|1) 
and 1 2)). This is easily done. 

So far there have been no approximations made. The 
essence of the dressed-state approximation is similar to 
the secular approximation in that we keep the internal 
dressed-states but allow the coherences to go to zero. In 
contrast to the secular approximation, we d o no t simply 
set pi2 = 0. Instead, we notice from Eq. (5_^) that, if 
we ignore the operator nature of the eigenenergies, then 
the equations for pi2 would have a term of the form 
— i (£'1 — E2) pi2- To first order, E1—E2 is approximately 
6. Thus pi2 will rotate very quickly such that only terms 
rotating at this very rapid pace will be able to contribute 
to its evolution. Thus only terms involving pi2 are kept 
in the equation for pi2 

P12 = —TB {sin 9 cos 9 p 12 cos 9 sin 9) 



^2 



3i2 + P12 sm fj 
-i [Ei{x,t)pi2 - Pl2E2ix,t)] 



(5.9) 

Knowing that the last term in Eq. ( |5.9| ) serves to force 
P12 to oscillate very rapidly and the other two terms force 
it to decay quickly, pi2 will average to zero. Thus we set 
P12 and P21 equal to zero in the population equations 
giving 



Pii 



P22 



TB (sin 9 cos 9pii sin 9 cos 9 

'9,pii} -i[Ei{x,t),pii 
cos 9p22 sin 9 cos 9 + cos' 



sin^ 6*^22 sin^ 0) 



4< 



(5.10) 



TB (sin 

r 

~2 



5ii cos 



{sin^9,p22] -i[E2{x,t),p22] 



(5.11) 
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I f we cou ld ignore the operator nature of the rates in Eqs. 
(5.9- 5.11 ), we would obtain the same rates as given by 
Dalibard and Cohen- Tannoudji in [ p^ . 

All that remains now is to design the master equation 
that forces pi2 to zero giving these population equations. 
The master equation which fulfills these requirements is 

p = -ICp-i[Eiix,t)\l){l\+E2ix,t)\2){2lp] 
+TBJ [cos sin {a^a - aa^)] p 
—TA [cos sin 6* (a^a — aa^)] p 
+r {BJ [cos^ ea]p~A[ cos^ 0a\ p) 
+r {BJ [sin^ ea^] p-A [sin^ 0a^] p) , (5.12) 

where the kinetic energy term has been restored. This 
form of the master equation is still very complicated re- 
membering the definitions of cos0 and sin0 from Eqs. 
(5.7). It is possible to simulate this master equation as 
it stands but the simulation would be very slow and in- 
efficient. This master equation can, however, be approx- 
imated further with almost no loss in accuracy. 

The definitions in Eqs. (5/7) can be approximated by 
remembering that we are working in the regime of S ^ SI. 
Thus, to leading order, cos 9 is simply 1, and sin 9 is n/26. 
Also any terms of order sin'* are extremely small and so 
are also ignored. If we propagate these approximations 
through our system, we find that the dressed-state |1) is 
very close to the bare excited state |e). Also, the dressed- 
state 1 2) is very near the bare ground state \g). Thus if 
we make the approximations 



|1> ^ \e) , |2) ^ Iff), 
then we can similarly approximate 



a ~ a, (7^ ~ a' . 



(5.13) 



(5.14) 



The last approximation is to expand the dressed-state 
eigenenergies in a Taylor series to second order giving 







4(5 


16(53 ' 






4(5 


16(53 J 



(5.15a) 



(5.15b) 



This leaves us with the final master equation for the 
dressed-state approximation as 



P = 



Px 

2m 



4(5 16(53 
^T{BJ[a]p~A[a] p) 

-T ( BJ 





P~ A 













(5.16) 



Here we have removed the Hq — Sa^ a Hamiltonian term 
by moving into an interaction picture. 

Again comparing this master equation to those already 
seen, the Hamiltonian terms here are exactly analogous 



to those der ived in the sophisticated adiabatic approxi- 
mation Eq. ( 3.13| ). Again we have the same spontaneous 
emission term as in Eq. (1.2b), but this time, the higher 
order correction term includes a spontaneous emission 
without changing the internal state of the atom, the op- 
posite of the case in the secular approximation. 

This would be a useful approximate master equation, 
especially if the atom was initially in the excited state. In 
our case however, we can simplify this equation further by 
noting that the jump terms keep the atom in the ground 
state. Once the excited state populations are reduced to 
zero, this master equation reduces to exactly the same 
master equation as derived for the more sophisticated 
adiabatic approximation Eq. (3.13). The dressed state 
master equation results would thus lie exactly on top of 
those of the more sophisticated adiabatic approximation 
and as such are not included in the simulations. 



VI. SIMULATION OF THE MASTER 
EQUATIONS 

Now that we have derived the forms of the master 
equations which we wish to compare, we need to set up 
a method of simulation for the different approximations 
and the full master equation. In this case, the numer- 
ical environment MATLAB turned out to be the most 
useful tool, combined with the Quantum Optics Toolbox 
produced by S.M. Tan |l^, |l^. The simulation is de- 
signed by converting states and operators to vectors and 
matrices. Making this conversion requires a number of 
different adaptations of the theoretical master equations. 

Firstly, we need to chose a form for the complex Rabi 
frequency operator, f2. The form of the Rabi frequency 
operator we use here is $7 = fimax sin kx, such that there 
is no time dependence (standing wave) and the opera- 
tor is Hermitian. Rewriting the sine function in terms of 
exponentials, e^**^^, we can examine the action of these 
exponentials. This suggests using the momentum repre- 
sentation, because the exponentials, e"^^^^ , simply repre- 
sent single momentum kicks to the atom of =plfifc, or =pl 
atomic unit of momentum. 

Having chosen to use the momentum representation 
to evaluate our master equation solutions, we need to 
address concerns with the conversion process. Firstly, 
the momentum range must be truncated. The loss of 
probability due to this truncation should be kept below 
some small level, say 10~^. Using this rule of thumb, it 
was found that a momentum Hilbert space ranging from 
—25hk up to +25hk was sufficient for the parameters we 
chose (discussed later). 

We still face a limitation problem in that the momen- 
tum Hilbert space is continuous. To simulate this on 
a computer, we need to discretize this space. We are 
able to discretize this space, again as a consequence of 
the choice of f2. The exponential component nature of 
51 provides for momentum kicks of exactly one atomic 
unit of momentum in the direction of propagation of the 
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laser light (the a;-dhection) . Thus the only momentum 
kicks that are not in a single unit of atomic momentum in 
the x-direction are due to spontaneous emission of pho- 
tons from the atom. This allows for a momentum kick 
of one atomic unit in any direction, which, when pro- 
jected onto the x-direction, allows for a random kick of 
anywhere between —1 and 1 atomic unit. However the 
relative infrequency of the spontaneous emission allows 
us to approximate this by a kick of —1, 0, or 1. 

This approximation requires us to convert the integral 
over the atomic dipole distribution to a discrete sum. 
Following Ref. ||, we let 




W{u)du — > V{u), (6.1) 



where the discrete function V{u) is obtained by making 
sure that it has the same normalization, mean momen- 
tum kick and mean squared momentum kick. This gives 
V{-1) = i, V{0) = f and V{1) = i. ^ 

We now need to define our matrix notation. The sim- 
ulations here will set h, k and m to 1 so that the kinetic 
energy operator can be represented as 




- ^ n^\n){nl (6.2) 



where |n) denotes a momentum eigenstate. The e~*'^^ 
operator in a momentum Hilbert space just acts as a 
raising operator given by 

25 

R= + (6-3) 

n=-25 

The e**^^ operator is just the lowering operator . For 
simplicity, we take the initial density matrix to be the 
zero momentum state |0), given by a 51 x 51 matrix of 
zeros with a 1 in the very centre. 

These are all the approximations required to allow the 
simulation of the adiabatic approximations. The full 
master equation and the secular and dressed-state ap- 
proximations require the internal state information as 
well. This is achieved by constructing the tensor product 
of the 51x51 momentum Hilbert with a 2 x 2 internal 
state Hilbert space. This is easily done using the quan- 
tum optics toolbox. 

The only problem remaining is to adopt a method of 
simulation. There are a number of different simulation 
methods available as part of the quantum optics tool- 
box as well as any number of methods available using 
regular ordinary differential equation techniques. The 
method we use here is a built-in function of the Quan- 
tum Optics Toolbox called odesolve. This allows op- 
tions to be specified for use on both smooth and stiff 
ODE problems. The simulations presented here using 
odesolve were checked using first, a hybrid Euler and 
matrix exponential method, and second, a modified mid- 
point method combined with Richardson extrapolation 



from Press et al ||T^ . The results all agreed well but the 
odesolve method was by far the fastest and easiest to 
use. 

We wish to simulate the experimental regime of 
Hensinger et al, where S l^max ^ F ~ Q'^^^/S. How- 
ever, the actual experimental parameters would be pro- 
hibitively time-consuming to simulate. This is both be- 
cause of the separation in time scales between the fastest 
(S) and slowest dynamics (Ffimax), and also because of 
the basis size required. The latter is determined by the 
fact that pf-^^^/2m must be larger than the effective po- 
tential drop, fl^^^/A5. If we were to use the parameters 
of the experiment of ^ then we would need a basis size of 
more than ±70hk. However, we can scale the parameters 
down and still preserve the regime of the experiments. 

As well as working in the regime S ^ Omax ^ F 
^max/*^' we require for the validity of the adiabatic ap- 
proximations that F be much larger than the oscilla- 
tion frequency near the potential minimum, Wqsc- In 
our scaled units (ft = fc = to = 1), the latter is of 
order •\/f^max/'^- On this basis, we have chosen para- 
maters of (5 = 10^, Omax = 2 X lO-^ and F = 200, leaving 
nl^^JS = 400, giving coosc = 20. For Rb as in Ref. [|, 
we have k = 27r/780nm and to — 1.42 x lO^^^kg, so in 
SI units, the frequency unit is hk'^/m — 4.821 x 10"* s~^. 
Note that the F we have chosen is not the true radiative 
decay rate for Rb. 

The scaled time unit can be given meaning by exam- 
ining the spontaneous emission rate. For each of the 
approximations, the fluorescence rate is (to leading or- 
der) TW^/AS'^. Here fi^ is a time-averaged effective value, 
which would be somewhat less than r^^^x- This means 
that after a time period of AS"^ /rV,f^^^, we would expect 
there to have somewhat under one spontaneous emission. 
This time is 2 scaled time units. There is, however, a lot 
of evolution occurring in that time period. To compare 
the approximations, we look in detail at a period from 
to 2 time units, and also look at the long time results at 
8 time units. 



VII. RESULTS OF THE SIMULATIONS 

In comparing the results of the simulations, we first 
compare the accuracy of the approximations, and then 
the resources required to perform the calculations. 



A. Accuracy 

To compare the accuracy of the simulations, we look 
at the momentum distribution of the atom as it evolves 
through time. The interesting components of this evo- 
lution are the probability to have momentum and 
the probability to have 1 atomic unit of momentum as 
the atom evolves through time. The other probabilities 
evolve similarly to one of these two. 
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FIG. 1: The probability to have zero momentum versus scaled 
time, a) Overall picture of the probability to have zero mo- 
mentum. Individual lines are not resolvable at full size, b) 
The heavy solid line is the full master equation, the dashed 
line represents the standard adiabatic approximation while 
the dot-dash line represents the secular approximation. The 
full master equation is indistinguishable from the more so- 
phisticated adiabatic and dressed state approximations even 
at this size, c) The thin solid line represents the dressed- 
state approximation and the sophisticated adiabatic approx- 
imation. 



FIG. 2: The probability to have 1 atomic unit of momentum 
versus scaled time, a) Overall picture of the probability to 
have 1 atomic unit of momentum. Individual lines are again 
not resolvable at full size, b) The heavy solid line is the full 
master equation, the dashed line represents the standard adi- 
abatic approximation while the dot-dash line represents the 
secular approximation. The full master equation is indistin- 
guishable from the more sophisticated adiabatic and dressed 
state approximations even at this size, c) The thin solid line 
represents the dressed-state approximation and the sophisti- 
cated adiabatic approximation. 



Firstly, we will look at the probability to have zero mo- 
mentum over a relatively short timescale. The approxi- 
mations are so close to the full master equation that, at 
full size, they are almost impossible to distinguish from 
the full master equation. Fig. ^a shows an overall picture 
of how the probability to have zero momentum evolves 
through time. 

Fig. |l]b zooms in on a section of the full size figure to 
illustrate the differences between the approximations. As 
we can see, the secular and standard adiabatic approxi- 
mation evolutions both significantly lead that of the full 
master equation in time. The dressed-state and sophis- 
ticated adiabatic approximations are very close to the 
full master equation solution even at this magnification. 
Fig |l]c zooms in even closer to try to distinguish the so- 
phisticated adiabatic approximation from the full master 
equation. 

This trend continues as we investigate the probabil- 
ity to have 1 atomic unit of momentum in Fig. |^. The 
inaccuracy of the secular and standard adiabatic approx- 
imations are again evident in Fig. The more sophis- 
ticated adiabatic approximation is very close to the full 
master equation evolution and is still indistinguishable 



at this magnification. 

Fig. ^ provides a means of comparing the sophisti- 
cated adiabatic approximation to the full master equa- 
tion evolution in detail. As we can see, the more sophis- 
ticated adiabatic approximation is again very close to the 
full master equation. 

The last visual comparison to make is to see how the 
evolution described by the approximations matches that 
of the full master equation after a very long time pe- 
riod. The time period that has elapsed in Fig. ^ is 8 
time units, after which, we would have expected there to 
be a number of spontaneous emissions. At this point in 
time, we compare the probability densities described by 
the approximations and the full master equation. It is 
evident from Fig. ^ that the standard adiabatic and the 
secular approximations are quite poor methods for simu- 
lating the full master equation evolution over a long time 
period. The main reason for this is the leading behaviour 
such as we see in Fig. |^a. Even at long times, the proba- 
bility to have zero momentum is still oscillating and the 
leading behaviour means that the standard adiabatic and 
secular approximations are not oscillating in phase with 
the full master equation. Thus, even though they follow 
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dressed-state and sophisticated adiabatic approximations 
are very close to the fuU master equation, even at this 
long time. 



B. Resources 

The main resource required to perform the simulations 
is time. Although each approximation has different mem- 
ory and processing power requirements, these needs are 
reasonably accurately reflected in the time each simula- 
tion takes to run. The times quoted in table ^ are the 
times required to calculate the set of results from to 8 
time units and are quoted in seconds. 



FIG. 3: Long time probability density. The solid line repre- 
sents the full master equation which is indistinguishable from 
the dressed-state and more sophisticated adiabatic approxi- 
mations at this magnification. The dashed line represents the 
standard adiabatic approximation while the dash-dot line is 
the evolution described by the secular approximation. It is 
worth noting here that due to our approximations, the prob- 
ability density is only defined for integer values of momentum. 
We have "connected the dots"to aid the eye. 




Momentum in atomic units 

FIG. 4: Long time probability density. The heavy solid line 
represents the full master equation evolution while the thin 
solid line represents the sophisticated adiabatic and dressed- 
state approximations. Again the secular and standard adi- 
abatic approximations are represented by the dot dash line 
and the dashed line respectively. 



TABLE I: Comparison of the times each simulation requires 



to run. 

Simulation Time (s) 

Full Master equation 617.7 

Standard Adiabatic Approximation 16.98 

More Sophisticated Adiabatic Approximation 18.31 

Secular Approximation 76.43 

Dressed State Approximation" 18.31 



"The time quoted here is idential for the more sophisticated adi- 
abatic approximation because the equations are identical for an 
atom initially in the ground state. 

As we can see from the results in table |, the adia- 
batic approximations are both clearly the fastest. The 
secular approximation is just over 4 times larger than 
the adiabatic approximations. This is not entirely un- 
expected. We would have expected at least a doubling 
in time by using the methods that included state infor- 
mation. The secular approximation is supposed to force 
the coherences to zero. Unfortunately, using our method 
to simulate this allows zero to be anywhere up to 10~^. 
This, although small, still has to be processed explaining 
the four-fold increase in time. One reason it is just over 4 
times the time required for the adiabatic approximations 
could be that after the evolution has been calculated, 
there is still a partial trace to be performed to obtain a 
solution of the same form as the adiabatic approxima- 
tions. We could have limited this time by simulating two 
coupled equations instead of a full master equation and 
then would probably have only doubled the time taken. 



roughly the same shape, they are not necessarily at the 
same point in the oscillation as the full master equation. 
What is even more striking though is that the secular 
approximation evolution seems to follow a slightly differ- 
ent trend for the probabilities to have larger momentum. 
The intriguing features of the secular approximation are 
discussed in Sec VIIl. 



To examine how close the sophisticated adiabatic and 
the dressed-state approximations are to the full master 
equation, Fig. || focuses on a smaller section to pro- 
vide a comparison, ft is evident from Fig. ^ that the 



VIII. DISCUSSION 

The standard and more sophisticated adiabatic ap- 
proximations take similar times to simulate, and are 
much faster than the full master equation. The only 
difference between them is an extra potential term. It 
turns out, though, that this Hamiltonian term is quite 
important in accurately describing the motion of the 
atom. While the standard approach evolves too quickly 
and leads that of the full master equation evolution, the 
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more sophisticated approach with the modified potential 
does not suffer this problem. The evolution described 
by this more sophisticated approach is very close to the 
full master equation, even at long times. Of course the 
dressed-state approximation offers the same accuracy as 
the sophisticated adiabatic approximation which may be 
useful if we wanted to simulate an initially excited atom. 

These successes contrast the results from the secu- 
lar approximation. As one can see from Figs. ^ and 
the secular approximation not only leads the full master 
equation solution but it also predicts a lower probabil- 
ity to have either zero or one atomic unit of momentum. 
This result is surprising because the secular approxima- 
tion master equation is quite similar to the others. In- 
vestigating this further, we find that the secular approxi- 
mation master equation simulation shows that the prob- 
ability to have 25 atomic units of momentum increases 
exponentially much faster than any of the other approx- 
imations. To analyse this in another manner, wc find 
that for the secular approximation Tr[p] falls off from 1 
exponentially much faster than any of the other approx- 
imations or the full master equation. 

If we consider the d ressed-state approximation mas- 
ter equation, Eq. ( ^.16 ), we notice that the Hamiltonian 
terms are, to leading order, the same as the final secular 
approximation master equation Eq. (4.3). The only ad- 
vantage the dressed state approximation has in terms of 
it's Hamiltonian is the presence of the higher order term. 
The leading order Lindbladian terms are also identical 
with only the higher order terms differing. The dressed- 
state master equation includes Lindbladian terms which 
give rise to a momentum kick to the atom (from the B su- 
peroperator) without a change in the internal state of the 
atom. The higher order secular approximation master 
equation Lindbladian term provides the opposite. Here 
there is an internal state change without a momentum 
kick to the atom. Thus it is quite unexpected that the 
secular approximation solution should differ so greatly 
from the dressed-state approximation solution as they 
are identical in the leading order terms. 

It is a fairly simple matter to investigate the effect of 
the differences between the two approximations by sim- 
ulating only parts of the master equations. If we only 
included the Hamiltonian terms from each master equa- 
tion, the difference would only be the higher order term 
in the dressed state master equation. This term has the 
same effect here as it does for the sophisticated adia- 
batic approximation, correcting the leading behaviour. 
It does not however account for the secular approxima- 
tion predicting lower probabilities as it does in Figs. ^ 
and y. Thus to investigate the difference between the 
higher order Lindbladian terms, we drop the B superop- 
erator from the dressed-state master equation and add it 
to the secular approximation master equation. This does 
not correct the problems with the secular approximation 
nor does it severely affect the dressed-state approxima- 
tion. This only leaves the inherent difference that the 
dressed-state master equation provides a correction term 
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FIG. 5: Ratio of Tr[peepl/2m] to TrfPpee] for the sophisti- 
cated adiabatic approximation. All other simulations look 
similar. 



without requiring a state change where the secular ap- 
proximation forces a state change in it's correction term. 
Thus we have to conclude that the best approximation to 
the full master equation involves a Lindbladian term that 
does not change the internal state of the atom. Lacking 
it, the secular approximation is the poorest. 

One other point to note in this discussion is one of the 
conditions on which the adiabatic elimination techniques 
are based. This is that Fpee be much larger than JCp^e- 
This will be satisfied if 



Tr 



'2m 



/rTr[pee] < 1. 



(8.1) 



For the parameter regime of this investigation, it is not 
obvi ous that this holds. Fig || shows how the ratio in 
Eq. ( ^.l|) evolves through time for the sophisticated adi- 
abatic approximation. Here we see that numerically, this 
ratio is around 0.1 which is of the same order as the other 
ratios (such as ftmax/S) which are required to be small 
for our approximations. 

Finally, we discuss the possibility of simulating with 
the true experimental parameters. This is difficult be- 
cause of the stiffness of the full master equation, and the 
basis size required for all methods. The latter problem 
can be avoided by using quantum trajectories ||l9|, |2^, pi] . 
Hensinger et al [g[ actually use quantum trajectory simu- 
lations based on the secular approximation master equa- 
tion. It is however, possible to convert any master equa- 
tion of the Lindblad form to a quantum trajectory simu- 
lation j2^. All of the approximate master equations we 
have developed here have been written in the Lindblad 
form and as such all of these could be simulated using 
quantum trajectories. 
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IX. CONCLUSION 

There are a number of theoretical models for the mo- 
tion of an atom as it interacts with a light field. This pa- 
per has investigated the possibility of using four different 
approximations as opposed to using the full master equa- 
tion to simulate an experimental system. Two have been 
widely used in the past and two have not. We have given 
a detailed explanation of the mathematical principles to 
perform each of these approximations on a fairly general 
system. We have also compared them numerically to to 
the true dynamics from the full master equation. 

In a regime of particular experimental interest, we have 
found that the most accurate results are obtained from 
two approaches that we have introduced here, a sophisti- 
cated adiabatic approach and a dressed-state approach. 
These give identical equations in the regime of interest. 



and in terms of resources, they are almost as fast to sim- 
ulate as the standard adiabatic approximation. This has 
been most used in the past, but deviates significantly 
from the true dynamics for long times. The other ap- 
proximation that has been used in the past, the secular 
approximation, is even poorer. On top of the failings of 
the standard adiabatic approximation, it takes longer to 
simulate and appears to produce anomalous momentum 
diffusion. 
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